# this function test whether the data follow zero inflated poisson distribution (ZIPD).

#---define the iteration fuction for MLE

lambdaML <- function(lm,m,n,n0){
res = m*(1-exp(-lm)) - lm*(1-n0/n)
return(res)
}


# ---ML fit to the zero inflated poisson distribution

fitzipd <- function(ct){
n0 = length(ct[ct==0])
n = length(ct)
m = mean(ct)
l_ml = uniroot(lambdaML,c(0.0001,10),m=m,n=n,n0=n0)$root
pi_ml = 1 - m/l_ml
return(c(l_ml, pi_ml))
}




#----generate expected frequency from model

zipd <- function(l,p,x){
    if(x == 0){res = p + (1-p)*exp(-l)}
    else {res = (1-p)*l**x*exp(-l)/factorial(x)}
    return(res)
}





#-----calcualte the chi2

chi2 <-function(ct){
    unqct = unique(ct)
    obs = numeric(length(unqct))
    pred = numeric(length(unqct))
    param = fitzipd(ct)
    for(i in 1:length(unqct)){
        obs[i] = length(ct[ct==unqct[i]])/length(ct)
        pred[i] = zipd(param[1],param[2],unqct[i])
    }
    chisq = sum((obs - pred)**2/pred)
    dof = length(unqct) - 2
    pvalue = 1- pchisq(chisq,dof)
#plot(obs,pred,xlab='observed frequency',ylab = 'predicted frequency',main=paste('Chi2',as.character(chisq))
    plot(obs,pred,xlab='observed frequency',ylab = 'predicted frequency')
    grid()
    abline(0,1)
    return (pvalue)
}


#------- test ------
ct = rpois(1000,1.)
ct = c(ct,numeric(100))


# ---data ----
 
b = read.table('argmb_09282014.csv',sep=',',header=T)
